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Abstract We model the evolution of the spin frequency’s second derivative iz and the braking 
index n of radio pulsars with simulations within the phenomenological model of their surface 
magnetic field evolution, which contains a long-term power-law decay modulated by short¬ 
term oscillations. For the pulsar PSR B0329H-54, a model with three oscillation components 
can reproduce its variation. We show that the “averaged” n is different from the instanta¬ 
neous n, and its oscillation magnitude decreases abruptly as the time span increases, due to 
the “averaging” effect. The simulated timing residuals agree with the main features of the 
reported data. Our model predicts that the averaged i/ of PSR B0329H-54 will start to decrease 
rapidly with newer data beyond those used in Hobbs et al.. We further perform Monte Carlo 
simulations for the distribution of the reported data in |f'| and |n| versus characteristic age 
Tc diagrams. It is found that the magnetic field oscillation model with decay index a = 0 
can reproduce the distributions quite well. Compared with magnetic field decay due to the 
ambipolar diffusion (a = 0.5) and the Hall cascade {a = 1.0), the model with no long term 
decay {a = 0) is clearly preferred for old pulsars by the p-values of the two-dimensional 
Kolmogorov-Smirnov test. 

Key words: stars: neutron — pulsars: individual (B0329H-54)—pulsars: general — magnetic 
fields 

1 INTRODUCTION 

The spin-down of radio pulsars is caused by emitting electromagnetic radiation and by accelerating particle 
winds. Traditionally, the evolution of their rotation frequencies v may be described by the braking law 
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where n is the braking index, and iY is a positive constant that depends on the magnetic dipole moment and 
the moment of inertia of the neutron star (NS). By differentiating Equation ( 1 ), one can obtain n in terms of 
several observables, n = vv/ii^. For the standard vacuum magnetic dipole radiation model with a constant 
magnetic field (i.e. K = 0), n = 3 (Manchester & Taylor 1977). Thus the frequency’s second derivative 
can be simply expressed as 

i> = 3v^ jv. (2) 

The model predicts i> > 0 and \ v\ should be very small. 

However, unexpectedly large values of i> were measured for several dozen pulsars thirty years ago 
(Gullahorn & Rankin 1978; Helfand et al. 1980; Manchester & Taylor 1977), and many of those pulsars 
surprisingly showed i/ < 0. Some authors suggested that the observed values of i/ could result from a noise- 
type fluctuation in the pulsar period (Helfand et al. 1980; Cordes 1980; Cordes & Helfand 1980). Based on 
the timing data of PSR B0329H-54, Demiahski & Proszyhski (1979) further proposed that a distant planet 
would influence i>, and the quasi-sinusoidal modulation in timing residuals might be caused by changes in 
pulse shape, precession of a magnetic dipole axis, or an orbiting planet. Baykal et al. (1999) investigated 
the stability of i) for pulsars PSR B0823H-26, PSR B1706-16, PSR B1749-28 and PSR B2021H-51 using 
their time-of-arrival (TOA) data extending to more than three decades. This analysis confirmed that the 
anomalous i> terms of these sources arise from red noise (timing residuals with low frequency structure), 
which may originate from external torques applied by the magnetosphere of a pulsar. 

The relationship between the low frequency structure in timing residuals and the fluctuations in pulsar 
spin parameters v, and v) is very interesting and important. We call both the residuals and the fluctuations 
the “timing noise” in the present work, since we will infer that they have the same origin. Timing noise for 
some pulsars has been studied for over four decades (e.g. Boynton et al. 1972; Groth 1975; Jones 1982; 
Cordes & Downs 1985; D’Alessandro et al. 1995; Kaspi, Chakrabarty & Steinberger 1999; Chukwude 
2003; Livingstone et al. 2005; Shannon &. Cordes 2010; Liu et al. 2011; Coles et al. 2011; Jones 2012). 
However, the origins of the timing noise are still controversial and there is still unmodelled physics to 
be understood. Boynton et al. (1972) suggested that the timing noise might arise from “random walk” 
processes. The random walk in v may be produced by small scale internal superfluid vortex unpinning 
(Alpar, Nandkumar (& Pines 1986; Cheng 1987a), or short time {t ^ 10 ms for the Crab pulsar) fluctuations 
in the size of the outer magnetosphere gap (Cheng 1987b). Stairs, Lyne & Schemar (2000) reported long 
time-scale, highly periodic and correlated variations in the pulse shape and the slow-down rate of the pulsar 
PSR B1828-11, which have generally been considered as evidence of free precession. The possibilities 
were also proposed that the quasi-periodic modulations in timing residuals could be caused by an orbiting 
asteroid belt (Cordes & Shannon 2008) or a fossil accretion disk (Qiao et al. 2003). 

Recently, Hobbs et al. (2010, hereafter H2010) carried out the most extensive study so far of long term 
timing irregularities of 366 pulsars. Besides ruling out some timing noise models in terms of observational 
imperfections, random walks, and planetary companions, some of their main conclusions were: (1) timing 
noise is widespread in pulsars and is inversely correlated with pulsar characteristic age Tc\ (2) significant 
periodicities are seen in the timing noise of a few pulsars, but quasi-periodic features are widely observed; 
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for a longer data span and the magnitude of |fi| for a shorter data span is much larger than that caused by 
the magnetic braking of the NS; and (4) the numbers of negative and positive v are almost equal in the 
sample, i.e. iVp « N^- Lyne et al. (2010) showed credible evidence that timing noise and v are corre¬ 
lated with changes in the pulse shapes, and are therefore linked to and caused by changes in the pulsar’s 
magnetosphere. 

Blandford &. Romani (1988) re-formulated the braking law of a pulsar as z> = —K{t)v^, which means 
that the standard magnetic dipole radiation is still responsible for the instantaneous spin-down of a pulsar, 
and vv/v^ ^ 3 does not indicate deviation from the dipole radiation model, but only means that K{t) is 
time dependent. Considering the magnetospheric origin of timing noise as inferred by Lyne et al. (2010), 
we assume that magnetic field evolution is responsible for the variation of K (f), i.e. K = in which 

A = is a constant, and R (~ 10® cm), I (~ 10^® g cm^), and 0 (~ 7r/2) are the radius, moment 

of inertia, and angle of magnetic inclination for the NS, respectively. We can rewrite Equation (2) as 

i) = 3iy^/v + 2iyB/B. (3) 

Since the numbers of negative and positive i> are almost equal, it should be the case that B quasi- 
symmetrically oscillates, and usually \2i'B/B\ ^ 3v^ jv. In addition, it can be noticed that pulsars with 
Tc < 10® yr always have i) > 3v^ jv (H2010); a reasonable interpretation is that their magnetic field decays 
(i.e. B < 0) dominate the field evolution for these “young” pulsars. 

Therefore, Zhang & Xie (2012a, hereafter Paper I) constructed a phenomenological model for the dipole 
magnetic field evolution of pulsars with a long-term decay modulated by short-term oscillations, 

B{t) = Bd{t){l + '^ki sm{cj)i + 271^)), (4) 

where t is the pulsar’s age, and ki, (f>i and Ti are the amplitude, phase and period of the i-th oscillating 
component, respectively. Bd{t) = in which Bq is the field strength at the age to, the index 

a = 0 means the field has no long-term decay, and it was found that a > 0.5 for young pulsars with 
Tc < 10® yr (see Paper I for details). Substituting Equation (4) into Equation (1), we get the differential 
equation describing the the spin frequency evolution of a pulsar as follows 

= -AB{tf. (5) 

In paper I, we showed that the distribution of i) and the inverse correlation of i) versus Tc could be 
explained well with analytic formulae derived from the phenomenological model. In Zhang &. Xie (2012b, 
hereafter Paper II), we also derived an analytical expression for the braking index (n) and pointed out that 
the instantaneous value of n of a pulsar is different from the “averaged” n obtained from the traditional 
phase-fitting method over a certain time span. However, this “averaging” effect was not included in our 
previous analytical studies; this work is focused on addressing this effect. 

This paper is organized as follows. In Section 2, we show that the timescales of magnetic field oscil¬ 
lations are tightly connected to the i> evolution and the quasi-periodic oscillations appearing in the timing 
residuals, and the reported data of PSR B0329H-54 are fitted. In Section 3, we perform Monte Carlo simula¬ 
tions on the pulsar distribution in the v — Tc and n — Tc diagram. Our results are summarized and discussed 
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2 MODELING THE v AND N EVOLUTION AND TIMING RESIDUALS OE PULSAR PSR 
B0329+54 

PSR B0329+54 is a bright (e.g. 1500 mJy at 400 MHz’), 0.71 s pulsar that had been suspected of possess¬ 
ing planetary-mass companions (Demiahski & Proszyhski 1979; Bailes, Lyne, & Shemar 1993; Shabanova 
1995). However the suspected companions have not been conhrmed and their existence is currently con¬ 
sidered doubtful (Cordes & Downs 1985; Konacki et al. 1999; H2010). Konacki et al. (1999) suggested 
that the observed ephemeral periodicities in the timing residuals for PSR B0329H-54 are intrinsic to this NS. 
H2010 believed that the timing residual has a form that is similar to other pulsars in their sample. They 
plotted |i>| obtained from the B0329H-54 data sets with various time spans (see Figure 12 in their paper). 
For data spanning ^ 10 yr, they measured a large and significant i), and found that the timing residuals take 
the form of a cubic polynomial. However, no cubic term was found for data spanning more than ^ 25 yr, 
and |i>| became significantly smaller. The reported periods of the timing residuals for PSR B0329H-54 are 
1100 days, 2370 days, and/or 16.8 years (Demiahski & Proszyhski 1979; Bailes, Lyne, &. Shemar 1993; 
Shabanova 1995). 

In order to model the i> evolution for pulsar PSR B0329H-54, we first obtain v{t) by integrating the 
spin-down law described by Equations (4) and (5), and then the phase 

$(/) = f ( 6 ) 

Jto 

Finally, these observable quantities, v, v and v can be obtained by fitting the phases to the third order of its 
Taylor expansion over a time span Tg, 

^{U) = $0 + - to) + ih(L - fo)^ -f - to)^. (7) 

We thus get v, v and v for Tg from fitting to Equation (7), with a certain time interval of phases ATint = 
10® s (interval between each TOA, i.e. ATint = fz+i — D- 

We adopt a goodness of fit parameter to show how well the model matches the data, i.e. = X) Y? = 
^ where the subscripts M and D refer to the model results and the reported data, and (Ji are the 

uncertainties in the reported data. In order to minimize we adopt the Simulated Annealing Algorithm 
(SAA) to reach a fast convergence and avoid being trapped in a local minimum, and we use a simulation 
based on the Markov chain Monte Carlo (MCMC) methods for the fitting to explore the whole parameter 
space. 

In the upper panel of Eigure 1, we show the reported and the best-fitting (simulated) results of |i>| 
for various Tg for PSR B0329H-54; the reported data are read from Eigure 12 of H2010. There are three 
oscillation components involved in the simulation, and a = 0 is taken from Equation (4). The obtained 
smallest value of is 9.1, with the number of degrees of freedom being 20, and all the best-fit parameters 
for the three oscillation components are listed in Table 1 . xf for each reported data point is also shown in 
the middle panel; in the bottom panel, we show the corresponding n with the same oscillation parameters 
obtained above. The braking index n = vv jiP' obtained directly from Equation (5) is called “instantaneous” 
n; similarly, what is obtained by fitting phase sets to Equation (7) is called “averaged” n. It can be seen that 
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Fig.l I i>\, Xi and n for PSR B0329+54. Top panel; reported and fitted \i>\. The values reported 
by H2010 are represented by large cross symbols (i> > 0) and large circles (i> < 0); the best¬ 
fitting results are represented by small cross symbols {i> > 0) and small circles (i> < 0); the 
horizontal dashed line represents i> = 3z>^/iz. Middle panel: the goodness of fit parameter Xi for 
the fit of I fil. It is shown that the three-component model fits the reported data quite well. Bottom 
panel; instantaneous (solid line) and averaged (crosses) values of n. The horizontal dotted line 
represents n = 0. 

the averaged n has the same variation trends as i/, since \Ai//i/\ ^ 10“® and |Ai>/i>| ~ 10“^ are tiny, 
compared to |Ai>/i>| ~ 1. The magnitude of the first period of the averaged n is close to the instantaneous 
one, but it decays significantly due to the “averaging” effect. 

Since both the fits with one and two oscillation components are not very good and are certainly rejected 
by the x^ test (e.g. the smallest of the two component simulation is 128), and x^ is not significantly 
reduced after setting the index a as a free parameter, we thus conclude that three oscillation components 
are necessary for the fitting the variation in \v\. 

We use the Pearson Correlation Coefficient p = Covariance(j(:,y) measure the covariance 

W Variance(X) Varian.ce(y) 

between the parameters, where X and Y are the parameters to be tested. We show the joint posterior 
probability distribution between each pair of parameters in Figure 2, with p labeled in each panel. For each 
of the three oscillation components, their phase (j) is completely coupled with their period T. All other 
parameters are well determined independently. 

The timing residuals, after subtraction of the pulsar’s v and v over 36.5 years for PSR B0329H-54, are 
also simulated with exactly the same model parameters used for modeling i>. In the simulation, the following 
steps are taken: 

(i) We get the model-predicted TOAs with ATjnt = 10® s using Equation (6) over 36.5 yr, with the 











6 


Y. Xie, S.N. Zhang, and J.Y. Liao. 


Table 1 Summary of all the best-fitting parameters. The first row lists the best-fitting values for 
all the parameters, and the second row lists their Icr errors. 


Parameters Ti 

T 2 

Ti 

ki 

fc2 

ks 

(j>i 

4>2 

4>3 

(yr) 

(yr) 

(yr) 

(lO-'*) 

(10-^) 

(10--^) 

(rad) 

(rad) 

(rad) 

Best-fitting 15.870206 

49.03738 

7.869207 

4.05 

1.89 

2.58 

0.496 

0.071 

1.937 

Ict error 5.1x10“® 

1.6 X 10“'‘ 

2.5 X 10“® 

0.11 

0.12 

0.21 

0.148 

0.36 

0.427 


(ii) By fitting the TOA set to 


$(f) = $0 + i^oit - to) + -vo{t - ^o)^ 


( 8 ) 


we get $ 0 . 1^0 and uq- 

(iii) Then the timing residuals after the subtraction of v and v can be obtained by 

rr '^{ti)-{<^o + Vo{U-to) + \vo{U-to)‘^) 

^ res ) — 

^0 


(9) 


In Figure 3, we plot the reported timing residuals (from fig. 3 of H2010) with crosses, and the simulated 
result of the model with three oscillation components with a solid line. Note that the simulated result is 
not the fit of the model to the reported timing residual. It is actually the application of the model with the 
parameters derived from the fitting of \i>\, i.e. the figure shows a comparison of the timing residuals of 
the model’s prediction with the reported data. The RMS of the reported residuals and the differences are 
0.0086 and 0.0048, respectively, i.e., nearly a factor of two reduction of timing noise in terms of RMS with 
the application of the three-component model. In order to show the effectiveness of the three-component 
model, we perform and F-test for the three-component model and the base model adopted in the TEMP02 
program (Hobbs et al. 2006). The F statistic is given by 

ix^ - xl)/idi - d2) 


F = 


Xl/d2 


27, 


( 10 ) 


where Xi and X 2 are Pearson values, i.e. = X) ■ 5 ^’ where Ri is the residual of the i-th point, and 
di = 133 and d 2 = 124 are the number of degrees of freedom for the base model and three-component 
model, respectively. Here we assume ai = ao, i.e, all data points have the same weight; this way, the result 
of the F-test is independent of the exact value of (Tq- F ^ 27 means that the probability to reject the three 
component model over the base model is less than 2.7 x 10“^®, and thus the significance of the three- 
component model over the base model is higher than lOcr. Our model implies that the timing residuals are 
also caused by the magnetic field oscillation, and the quasi-periodic structures in timing residuals have the 
same origin (which is determined by Equation (5)) as those in i>, i> and ly variations. On the other hand, the 
fit is worse as the time span increases up to ~ 30 years, which may be mainly due to the additional noise 
components not included in i/ variations. 

The model includes an oscillation component with a period of ~ 49 years, however, it is hard to test 
directly from the power spectrum of its timing residuals, since the period is longer than the observational 
data span. However, there are still some features demonstrating its existence. Eor instance, the observed 
data were reported about four years ago, and the model predicts that F of the pulsar is now experiencing 
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Fig. 2 Correlation of the nine fitting parameters of \i/\. Each panel shows the joint posterior 
probability distribution between a pair of parameters, with correlation coefficient p labeled in 
it. The oscillation periods Ti, T 2 , and T 3 are plotted in terms of the differences from their mean 
values Tim = 5.00482804 x 10® s, Tam = 1.546443564 x 10® s, and Tgm = 2.48163290 x 10® s, 
respectively. 


data. The test could also be conducted by applying the model to a larger set of pulsars, which have short 
oscillation periods (shorter than the observed time span), and relatively large oscillation amplitudes (so that 
the swing behavior of could emerge; the exact criteria of k depend on v, v and T). 

From Equation (5), we can obtain an analytic approximation from the one oscillation component model 
(in Paper I) for i) 

O.jrt 

i>~- 2 i>(a/t + /cos(—+ ()))), ( 11 ) 

where / = 27 rfc/T represents the magnitude of the oscillation term. Thus, both parameters k and T are 
important. For a = 0, Equation (11) can be simply rewritten as 

27rf 

i) ~-2vf cos{—+ (j)), ( 12 ) 


One can see that the model predicts an oscillation behavior of i>, which implies that one may get either a 
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Fig. 3 Timing residuals of PSR B0329+54. The reported timing residual, after subtraction v and 
V of the pulsar over the 36.5 years, is represented by crosses; the predicted residuals modeled by 
three oscillation components are represented by the solid line. The model parameters are identical 
to those for the v simulation, as shown in Figure 1 . 

From Table 1 , we know that /i « /a ^ for PSR0329+54. Therefore the second component is 
less important in contributing to v, according to Equation (12). It is possible that the third component is the 
higher harmonic of the first component, since Ti Ri 2 T 3 . Thus it is likely that the first oscillation component 
dominates the timing behavior of the pulsar. As a matter of fact, there is almost always one dominant peak 
in the power spectrum of the timing residuals of most radio pulsars (H2010), i.e., one dominant oscillation 
component associated with their magnetic evolution. 

A major prediction of this model with three oscillation components is that the averaged i> will start to 
decrease rapidly with additional data that extend just a few years beyond the span that was used in H2010, as 
shown by the black crosses in the upper panel of Figure 1. As the data are already available to the observers, 
we suggest that this prediction can be used to confirm or deny our model. 

3 SIMULATING THE DISTRIBUTIONS OF i) AND N AND THEIR CORRELATIONS WITH re 

In this section, based on our phenomenological model, we use the Monte Carlo method to simulate the 
distributions of i/ and n, and their correlation with Tc. The “averaging” effects are naturally included in the 
simulations. For simplicity, in the following simulation we assume that there is one dominant oscillation 
component, which mainly determines the variations of i> and the timing residuals, as discussed above. If the 
one oscillation component model is rejected by the reported data, then a multiple component model shall 
be presented. This is not in conflict with the above three components fit, since fitting to the distributions 
requires much less detailed information about variations in i> for individual pulsars. 

We assume that the sample of phases (p of the field oscillation follows a uniform random distribution 
in the range from —tt to tt. Randomly Drawing a data set {i/, i>, Tg} from the reported sample space, i.e. 
from table 1 of H2010, calculating a corresponding start time to, and assuming some certain values for k 
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for TOAs is also assumed to be a constant, i.e. ATi„t = 10® s. Then the “averaged” values of v, v and v 
can be obtained by fitting to Equation (7). Hence one has its averaged |i>|, |n| and Repeat this 

procedure TV times, we will have N data points in \ v\-Tc and \n\-Tc diagram. 

3.1 Effects of Oscillation Period and Amplitude 

Analysis of a large sample of pulsar timing noise (H2010) showed that the oscillation periods are usually 
on the order of about 10 years. However, the structures seem to vary with data span, and as more data 
are collected, more quasi-periodic features are observed. In this subsection, we investigate the ranges of 
variation of k for a series of T. 

We show the measured |i>| and |n| versus Tc for 341 normal radio pulsars with 10® yr in Figure 
4, in which 184 pulsars with positive i> and n are plotted in the left panels, and the other 157 pulsars with 
negative values in the right panels. The simulated results, for the case of T =10 yr, with /3 = 4.6 or 2.1 
are also shown in the panels, in which j3 is defined hy k = 10“^ for convenience. One can see that the 
envelopes of ^ = 4.6 and 2.1 lie around the lower and upper boundaries of the reported data, respectively. 
This gives a natural constraint for k. Meanwhile, in the simulation the number of data points with > 0 
should be roughly equal to the number of < 0, i.e. Np/N^ ~ 1. In Table 2 we summarize the ranges of 
variation of (3 for different T. Physically, /3 < 0 is unacceptable, thus it should be T < 10® yr. However, our 
model fails to give a tight constraint on T. Note that Tc in the figure is the characteristic age of the pulsars. 
However, Tg in Figure 1 is the time span of the observation. Thus, the positive correlation between Tc and 
n in Figure 4 is not in conflict with Figure 1 which shows n decaying with Tg. 

Table 2 Ranges of variation of f3 for different T. /3„iin and ^max are the minimum and maximum 
values of P, respectively. 


T (years) 

10 

100 

1000 

10"* 

10® 

(/3min) /3max) 

(2.1, 4.6) 

(1.9, 4.5) 

(0.8, 3.4) 

(0.3, 2.3) 

(-0.5, 1.6) 


3.2 Two-dimensional Kolmogorov-Smirnov Test 

In this subsection, we perform a two-dimensional Kolmogorov-Smirnov (2DKS) test to reexamine the con¬ 
sistency of the distributions of the reported data and the simulated data, using the KS2D package . If the 
returned p-value is greater than 0.2, then it is a sign that we can treat them as drawn from the same distri¬ 
bution. 

We regard T as a random number and allow it to vary from 1 to 100 years to account for the diversity 
of periodicities observed in the population. Fet /3niin and /3max vary from 1.5 to 2.5 and from 4.0 to 5.0, 
respectively. It is found that (3 varying from 2.1 to 4.5 gives the highest p-value, as shown in the upper four 
panels of Figure 5. The returned probabilities are also labeled in each panel. Since the p-values indicate 
that the two samples are highly consistent, we thus conclude that the one oscillation component model with 
a = 0 is good enough to reproduce the liil-Tc and \n\-Tc distributions. 
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Fig. 4 Simulations of the \ v\-Tc distribution (upper panels) and \n\-Tc distribution (bottom panels) 
for r = 10 years. 


The simulated distributions with different a are also examined with a 2DKS test. We show the simulated 
distributions with a = 0.5 and 1.0 in the middle and bottom four panels of Figure 5, respectively. To 
describe the evolution of the pulsar magnetic field, three routes are generally proposed (see e.g Goldreich 
& Reisenegger 1992), i.e. ohmic dissipation. Hall effect, and ambipolar diffusion. Power law decay with 
a = 0.5 and 1.0 are produced by ambipolar diffusion and Hall effect, respectively (Paper I). Here we do 
not include ohmic dissipation since it is not important for pulsars with r > 10^( g (Gumming et 

al. 2004). For both cases of a = 0.5 and 1.0, the p-values of the simulated data for F < 0 are much lower 
than 0.2, and thus are rejected by the test. In fact, one can see that for Tc > 10® yr there is a crowded area of 
data points along the lower boundary for i> > 0, and the data points are scarce around the lower boundary 
for i> < 0. This is mainly caused by the long-term magnetic field decay, i.e. the decay term —2z>Q:/f > 0 
dominates the oscillation term —2z>/cos(^ + (j>) in Equation 11 for some cases. However, there is no 
such crowded area or scarce area in the reported data, which clearly indicates that the model with a = 0 is 
preferred for pulsars with Tc > 10® yr. 


4 SUMMARY AND DISCUSSION 

In this work we first modeled the i> and n evolutions and applied the obtained model parameters to simulate 
the timing residuals for the individual pulsar PSR B0329+54. Using a Monte Carlo simulation method, we 
simulated the distributions of pulsars in the |F| — Tc and |n| — Tc diagrams, and compared the simulation 
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Fig. 5 Simulations of \iy\-Tc and \n\-Tc distributions for T varying randomly from 1 to 100 years. 
The cases of a = 0, a = 0.5 and a = 1 are shown in the upper four panels, middle four panels 
and bottom four panels, respectively. 


1. We modeled the i) evolution of pulsar PSR B0329+54 with a phenomenological model that incorporates 
the evolution of B, which contains three oscillation components (upper panels of Figure 1). The model 
can reproduce the |i>| variation quite well, including the swings between i) > 0 and F < 0. This model 
predicts that the averaged i> of PSR B0329+54 will start to decrease rapidly with newer data beyond 
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2. We showed that the “averaged” values of n are different from the instantaneous values (bottom panels of 
Figure 1), and the oscillation abruptly decays after the first period due to the “averaging” effect. Using 
these parameters obtained from modeling the evolution of the averaged i>, we simulated the timing 
residuals of the pulsar (Figure 3), which agrees with the reported residuals (H2010) well. 

3. We performed Monte Carlo simulations for the distribution of |j>| and |n| in the \i>\ — Tc and \n\ — 
diagrams, respectively. The simulated results for different modes of long-term decay of the magnetic 
field (i.e. a = 0, a = 0.5 and 1.0 in Figure 5) are tested by the 2DKS. It is found that the reported 
distributions can be well reproduced with the one-oscillation-component model with a = 0 for pulsars 
with Tc > 10® yr. 

Pons et al. (2012) proposed a similar model of magnetic field oscillations with a timescale of (10® — 
IqS) io__G yj. magnitude 6B/B ~ 10“®, and obtained pulsar evolutionary tracks in the P — P diagram. 
Lyne et al. (2010) showed credible evidence that timing residuals and F are connected with changes in the 
pulse width. Therefore, timing residuals are more likely caused by changes in a pulsar’s magnetosphere 
with periods of about 1 — 100 yr. In Xie & Zhang (2013), we suggested that perturbations from Hall waves 
in the dipole magnetic field associated with NS crusts are probably responsible for the observed quasi- 
periodic oscillations in the timing data as well as changes in the pulse width, which may provide a physical 
explanation for the present model. 

We therefore conclude that magnetic field oscillations dominate the long term spin-down behaviors 
of old NSs, for which the long-term field decay is not important, in contrast to younger NSs with Tc < 
10® yr. The fact that only one oscillation component is required to reproduce the observed |F| — Tc and 
|n| — Tc distributions suggests that there is one dominant oscillation component for most NSs, and thus 
does not conflict with the fact the multiple oscillation components are also often observed in some pulsars. 
In fact, for some pulsars, the structures seen in the timing noise vary with data span and more quasi- 
periodic features are observed for longer data span (H2010). Admittedly, our current model cannot predict 
the number, amplitudes and periods of oscillation modes. However, our model can adequately describe the 
acquired timing data with a small number of oscillation modes, as shown in section 2, which represents a 
first step towards understanding of the magnetic field oscillations of NSs. As such, our understanding of the 
oscillation modes will be improved as more quasi-periodic features are revealed with longer observations 
in the future. In addition, our model can also describe the distributions of F and n reasonably well. As far 
as we are aware of, our work is the first one in which the distribution of F is used to test the long-term 
magnetic field evolution of NSs, which is independent from tests based on the traditional — F diagram. 
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